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A  PARAMETRIC  REPRESENTATION  OF  A  THUNDERSTORM  CELL 

1.  INTRODUCTION 

Rapid  advances  in  computer  technology  within  the  last  decade  have  led  to  dramatic 
increases  in  the  use  of  modeling  and  simulations  within  all  branches  of  the  Department  of 
Defense.  Models  and  simulations  now  are  utilized  in  many  functional  areas,  including  training 
and  military  operations,  research  and  development,  test  and  evaluation,  and  cost  analysis. 
Simulations  which  deal  with  atmospheric  effects  or  require  environmental  data  may  utilize  feature 
models,  which  provide  simple  parametric  representations  of  specific  atmospheric  features.  A 
feature  model  may  be  superimposed  on  an  initial  data  field  to  increase  atmospheric  fidelity  or 
provide  time  variability,  or  may  be  used  as  a  unique  environmental  source  for  the  simulation. 

This  report  presents  a  feature  model  of  a  thunderstorm  cell,  the  basic  constituent  of  the 
common  multi-celled  thunderstorm.  Unlike  numerical  convective  models,  which  describe  the 
development  and  structure  of  a  thunderstorm  using  mathematical  equations  to  represent  its 
dynamic  and  microphysics  processes,  the  feature  model  described  herein  (TSTORM)  only 
provides  a  simplistic  representation  of  thunderstorm  structure  based  on  equations  which  combine 
elementary  algebraic  and  transcendental  functions.  While  this  parametric  representation  of  a 
thunderstorm  is  not  as  physically  consistent  as  a  numerical  model  and  does  not  provide 
developmental  and  structural  aspects  possible  with  the  more  complex  numerical  model  (e.g.,  a 
storm's  life  cycle,  water  phase  changes,  heat  transfers),  it  does  provide  ease  of  use,  flexability  and 
implementation  capability  on  the  simulator  computer  system. 

As  a  data  source  to  a  host  simulation,  the  TSTORM  algorithm  is  capable  of  providing 
a  three  dimensional  (3-D),  steady-state  depiction  of  several  meteorological  fields  within  a 
"typical"  mature  thunderstorm  cell;  these  fields  include  the  vertical  velocity,  the  total  water 
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content  (liquid  plus  solid),  the  horizontal  wind  direction  and  speed,  and  the  temperature  and 
pressure  anomalies  relative  to  the  non-storm  environment.  Electrical  activity  and  lightning 
within  the  thunderstorm  cell  are  not  modeled.  A  dropsize  distribution  based  on  an  empirical  rela¬ 
tion  from  Marshall  and  Palmer  (1948)  is  available  at  user-specified  storm  locations.  Variations 
in  storm  intensity  are  provided  in  TSTORM  through  adjustments  in  a  few  initialization 
parameters.  If  superimposed  on  a  coarser  data  field,  the  time  history  of  a  thunderstorm  cell  may 
be  simulated  by  propagating  TSTORM  environmental  fields  past  a  selected  location. 

2.  TSTORM  DESCRIPTION 

2.1  Design  Framework 

The  TSTORM  algorithm  computes  meteorological  fields  within  a  3-D  (x,y,z)  framework 
(see  Figure  la).  Here,  both  the  x-  and  y-axes  extend  over  the  interxal  [-1,1],  while  the  vertical 
axis  extends  from  the  ground  to  the  cloud  top  (H).  In  TSTORM,  the  y-  and  z-coordinates 
represent  distance;  the  x-coordinate  can  either  represent  distance  or  time.  If  the  x-axis 
corresponds  to  distance,  then  TSTORM  provides  a  3-D  steady-.siale  depiction  of  a  thunderstorm 
cell;  for  this  representation,  the  points  A  and  B  of  Figure  la  would  be  surface  locations  on  the 
storm's  perimeter  while  point  C  would  be  at  the  storm's  (coordiiiate)  center.  On  the  other  hand, 
if  the  x-coordinate  represents  time,  a  (y,z)  vertical  slice  at  any  selected  value  of  x  on  the  interval 
[-1,1]  provides  an  instantaneous  representation  of  the  simulated  thunderstorm  cell  structure  at  a 
specific  time  during  passage  of  the  cell  at  a  station.  The  designation  of  the  x-axis  as  time  is 
possible  since  TSTORM  parameterizations  are  developed  on  the  assumption  that  the  storm's 
translational  motion  (V-j)  is  along  the  x-axis.  In  this  case,  the  points  A  and  B  of  Figure  la 
correspond  to  pre-  and  post-storm  locations,  respectively,  while  point  C  corresponds  to  the  time 
when  the  storm  is  directly  overhead. 
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The  TSTORM  meteorological  fields  are  computed  on  a  3-D  grid  (see  Figure  lb).  Here, 
both  the  X-  and  y-dimensions  are  set  to  21,  corresponding  to  a  grid  spacing  of  0.1  units.  The 
vertical  coordinate  NLVL  (the  number  of  levels)  is  of  variable  dimension  depending  on  what 
value  of  H  is  selected  (NLVL  =  2*H  +  1);  the  TSTORM  algorithm  presets  the  vertical  grid 
spacing  to  0.5  km. 

2.2  Inputs,  Outputs,  and  Limits 

With  reference  to  the  computer  program  listing  provided  in  the  Appendix,  inputs  and 

outputs  to  the  TSTORM  subroutine  subprogram  are  executed  within  the  calling  program  (i.e.,  the 

demonstration  driver  program)  by  means  of  the  following  statement: 

CALL  TSTORM(NLVL,H,WMAX,VMAX,WDSFC,VEER.NLOC,XDSD,YDSD,ZDSD, 
W,TA,PP,WC,WS,WD,DD,DSD,RR) 

The  input  parameters  (and  their  permissible  values)  which  are  passed  to  the  TSTORM  subroutine 
and  which  need  to  be  specified  by  the  operator  are: 

1)  Number  of  levels  (NLVL)  -  13  to  31 

2)  Thunderstorm  height  (H)  -  [6.0,6.5,...14.5,15.0]  km 

3)  Storm  maximum  vertical  velocity  (WMAX)  -  6  to  15  m/s 

4)  Wind  speed  at  level  H  (VMAX)  -  10  to  30  m/s 

5)  Surface  wind  direction  (WDSFC)  -  0  to  360  deg 

6)  Wind  veer,  surface  to  level  H  (VEER)  -  +  360  deg 

7)  Number  of  sites  for  dropsize  distribution  (NLOC)  -  >  1 

8)  Coordinates  (x,y,z)  for  each  dropsize  distribution  site 

XDSD(NLOQ  -  [-1.0,-0.9,...0.9,1.03  units 
YDSD(NLOQ  -  [-1.0,-0.9,...0.9,1.0]  units 
ZDSD(NLOC)  -  [0.0,0.5,...4.5,5.0]  km 
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The  parameters  H  and  WMAX  are  very  important  since  they  control  the  shape  and 
intensity  of  the  thunderstorm  simulation.  Generally,  the  larger  the  values  of  H  (storm  depth)  and 
WMAX  (maximum  storm  vertical  velocity),  the  more  severe  the  thunderstorm.  The  range  in 
permissible  values  for  H  and  WMAX  allows  the  operator  to  simulate  thunderstorm  cells  of 
varying  intensity,  from  light  to  moderate,  heavy  ic  severe.  Note  that  TSTORM  requires  H  to  be 
specified  as  a  whole  number,  or  a  whole  number  plus  0.5. 

The  fourth  through  sixth  input  parameters  specify  the  background  or  environmental 
(outside  storm)  wind  characteristics.  The  parameter  VMAX  is  important  since  it  controls  the 
entire  environmental  wind  speed  profile,  from  H  to  the  surface;  the  higher  (lower)  the  value  of 
VMAX,  the  higher  (lower)  the  wind  speed  throughout  the  vertical  column.  Meteorologically, 
wind  veer  is  defined  as  positive  (negative)  for  wind  which  rotates  clockwise  (counterclockwise) 
with  height.  For  example,  a  surface  wind  from  the  east  (90  deg)  which  veers  a  positive  135  deg 
will  be  a  southwesterly  wind  of  225  deg  at  top  level. 

Input  parameters  seven  and  eight  allow  the  user  to  specify  any  number  of  thunderstorm 
locations  for  dropsize  distribution  (DSD)  data.  Since  the  TSTORM  algorithm  does  not  provide 
for  spatial  interpolation,  these  DSD  locations  must  be  input  in  Cartesian  (x,y,z)  coordinates  which 
correspond  exactly  to  TSTORM  grid  locations.  For  example,  a  distribution  location  specified  as 
(-0.32,0.18,1.2)  would  be  unacceptable,  but  (-0.3,0.2,1.0)  would  be  valid.  In  TSTORM,  dropsize 
distribution  is  computed  using  an  empirical  relation  based  on  "rain"  drop  measurements.  Thus, 
in  order  to  increase  the  likelihood  that  dropsize  distributions  will  be  specified  for  liquid  (and  not 
solid)  water  drops,  TSTORM  limits  DSD  locations  to  below  5  km. 
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After  TSTORM  subroutine  execution,  the  following  six  gridded  meteorological  fields 
(dimensioned  (21,21,NLVL))  are  returned  to  the  demonstation  driver  program  and  are  available 
to  the  operator; 

1)  Vertical  velocity  W  (m/s) 

2)  Temperature  anomaly  (C) 

3)  Pressure  perturbation  Pp  (mb) 

4)  Total  water  content  WC  (g  m'^  ) 

5)  Wind  speed  WS  (m/s) 

6)  Wind  direction  WD  (deg) 

Also  returned  to  the  calling  program  are  the  drop  diameter  array  DD  (containing  21  diameters 
D  from  0  to  5  mm),  the  distribution  of  the  number  of  rain  drops  versus  diameter,  per  millimeter 
diameter  interval  and  per  cubic  meter  of  air  (DSD),  and  the  rainfall  rate  RR  (in  mm/hr).  Both 
DSD  and  RR  are  specified  for  each  dropsize  distribution  location  requested  by  the  operator. 

The  TSTORM  temperature  and  pressure  fields  (T,  and  Pp )  are  both  specified  as  in-storm 
anomalies  relative  to  the  background  (i.e.,  outside  storm)  environment.  In  order  to  define  an 
actual  temperature  (pressure)  structure  within  the  thunderstorm,  the  temperature  anomaly  (pressure 
perturbation)  field  would  be  added  to  the  environmental  temperature  (pressure).  Although 
TSTORM  does  not  provide  background  temperature  or  pressure,  they  could  be  provided  by  the 
simulator  in  the  form  of  actual  temperature  or  press"re  vertical  profiles.  Alternatively,  a 
background  temperature  profile  could  be  specified  by  a  temperature  at  a  particular  height  (usually 
the  surface)  and  an  assumed  environmental  lapse  rate  (typically  5  to  6.5  C/km).  The  hypsometric 
formula,  which  relates  height,  temperature  and  pressure,  could  then  be  used  to  determine  the 
background  vertical  pressure  profile. 
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In  TSTORM,  both  initialization  and  output  parameters  are  generally  specified  in  the  metric 
system  of  units;  the  algorithm  has  no  capability  for  conversion  from  metnc  to  other  systems  of 
units  (e.g.,  English).  Additionally,  the  TSTORM  algorithm  docs  not  provide  graphical  output 
displays;  these  are  both  user-  and  machine-dependent  and  are  left  to  the  discretion  of  the 
simulator.  The  demonstration  program  listed  in  the  Appendix  provides  the  tabular  data  output 
of  Figures  2  through  7,  and  14,  of  this  report. 

2.3  Parameterizations 

The  TSTORM  algorithm  produces  a  3-D  structure  of  significant  meteorological  parameters 
within  a  thunderstorm  cell  based  on  combinations  of  algebraic  and  transcendental  functions.  The 
vertical  velocity,  the  temperature  anomaly,  the  pressure  perturbation  and  the  total  water  content 
are  all  initially  calculated  at  grid  points  of  the  (x,z)  plane  above  the  x-axis  (i.e.,  at  y  =  0).  then 
expanded  horizontally  in  the  y-direction  to  obtain  a  full  3-D  representation.  The  horizontal  w  ind 
speed  and  direction  are  computed  on  the  (x,y)  plane  one  level  at  a  time,  starting  at  the  surface. 
The  dropsize  distribution  (DSD)  is  given  by  the  empirical  relation  of  Marshall  and  Palmer  (1948). 
Vertical  Velocity  W 

The  3-D  vertical  velocity  field  created  by  TSTORM  is  dominated  by  two  mid-level  cores  - 
an  updraft  core  in  the  forward  part  of  the  storm,  and  a  downdraft  core  w  ithin  the  rear  sector.  For 
given  values  of  x  and  z  at  y  =  0,  the  vertical  velocity  W(x,0,z)  is  obtained  as  the  product  of  the 
following  two  equations; 

Wl^e  *  sm(7i  *XDPL)  ^ 

and 


W2  = 


AWMAX 


H 


Z-- 


H 


1  + 


H) 


(2) 
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where 


kDPL^x 


(3) 


and  H2  H/2. 

The  curve  defined  hy  eq.  1  is  distinguished  by  its  sinusoidal  shape.  At  z  =  H2,  where 
XDPL  =  X,  eq.  1  provides  positive  numerical  values  (i.e.,  upward  motion)  for  all  x  >  0,  and 
negative  values  (for  downward  motion)  for  x  <  0.  The  parameter  XDPL  (eq.  3)  modifies  values 
of  X  to  provide  a  negative  vertical  tilt  (most  pronounced  at  small  and  large  values  of  z)  to  the 
eq.  1  curve.  The  parabolic  function  in  eq.  2  was  used  by  Kessler  (1969)  to  represent  the 
variation  of  vertical  velocity  with  height;  this  function  yields  a  maximum  value  at  mid-level  (H  2) 
and  values  of  zero  at  both  the  surface  and  top.  When  this  function  is  multiplied  by  eq.  1,  the 
desired  two  lobe,  mid-level  structure  for  W  is  obtained.  The  second  factor  of  eq.  2  is  used  to 
adjust  the  positive  W  lobe  slightly  upward,  and  the  negative  W  lobe,  slightly  downward. 

The  full  3-D  representation  of  the  vertical  velocity  is  given  by 

w:w)  =  *  YFAC 

where 

YFAC-{\-\y\f  for  >  =  [-1,1]  (5) 

The  function  YFAC  provides  exponential  decreases  in  vertical  velocity  values  away  from  the 
x-axis,  for  both  positive  and  negative  y  values. 
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Temperature  Anomaly 

The  TSTORM  3-D  depiction  of  the  temperature  anomaly  field  is  dominated  by  two 
opposing  cores  -  a  warm  core  located  aloft  below  the  updraft  maximum,  and  a  stronger  cold  core 
which  expands  and  strengthens  toward  the  surface.  For  given  values  of  x  and  y  at  y  =  0.  the 
temperature  anomaly  T^  (x.O.z)  is  determined  by  the  product  of  eq.  1  and  the  following  two 
equations: 

T1  *  (0.3///2)  ^  {H-z? 

and 

T2  =  1  *|x|)  /  /  24)} 


Additionally,  for  a  (x,  0.  z)  grid  point  located  within  the  rear  Ic.wer  sector  of  the  storm  (i.e., 
x  <  0  and  z  <  H2).  the  above  resultant  product  is  multiplied  by 


T3 


*  X  * 


4.32 


(8) 


Identical  to  the  vertical  velocity  parameterization,  the  temperature  anomaly  computation  starts 
with  the  sinusoidal  curve  of  negative  vertical  tilt  defined  by  eq.  1.  An  important  component  of 
the  Tg  computation  is  eq.  6,  w'hich  provides  large  numerical  values  to  the  computation  near  the 
surface  near  x  =  0,  and  zero  values  (at  all  x)  at  the  top  (z  =  H).  Equation  7  is  included  in  the 
Tj  computation  to  adjust  unreasonably  low  numerical  values  at  large  z  and  at  values  of  x  well 
removed  from  the  z-axis.  Another  adjustment  factor,  eq.  8,  is  included  in  the  computation  to 
insure  that  negative  temperature  anomaly  values  of  significant  magnitude  occur  near  the  surface 
for  all  X  <  0. 
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Analogous  to  the  vertical  velocity,  the  3-D  representation  of  the  temperature  anomaly  is 
given  as 

=  (9) 

where  YFAC  is  the  exponential  function  defined  by  eq.  5. 

Pressure  perturbation  Pp 

In  a  dynamic  sense,  the  pressure  perturbation  field  within  a  thunderstorm  is  linked  to  both 
the  vertical  velocity  and  temperature  anomaly  fields.  In  TSTORM,  these  two  fields  are  used  to 
create  a  pressure  perturbation  field  characterized  by  two  opposing  lobes  at  mid-storm  level,  and 
a  positive  anomaly  at  the  surface,  below  the  downdraft.  For  given  values  of  x  and  z  at  y  =  0, 
the  pressure  perturbation  is  computed  as 

P.ixM  =  0.1  *  -  r^u,o^) }  (10) 

for  z  <  HB  or  z  >  HT  ,  and 

for  z  >  HB  and  z  <  HT  , 
where  HB  =  H/12  and  HT  =  H  -  HB. 

Note  that  the  pressure  perturbation  disappears  at  cloud  top  height  H,  since  both  W  and  T,  are 
zero.  Analogous  to  the  vertical  velocity  and  temperature  anomaly,  the  full  3-D  representation  of 
the  pressure  perturbation  is  given  as 

/’/W)  =  Ppix,0^)  *  YFAC  (12) 
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Total  Water  Content  \VC 


The  parameter  WC,  which  includes  liquid  and  solid  water  hut  not  water  vapor,  provides 
the  cloud  form  as  w’ell  as  the  location  and  intensity  of  precipitation  for  the  TSTORM  simulation. 
The  main  features  of  the  WC  parameterization  include  a  mid-level  core  of  high  water  content, 
a  pronounced  near-surface  "rain"  shaft  near  the  storm's  center,  and  a  skewness  in  water  content 
at  mid-  to  upper  levels  in  the  direction  of  storm  motion.  Mathematically,  the  water  content 
WC(x,0,z)  is  initially  specified  by  the  addition  of  the  following  two  equations: 

WCl  =|7;(ar.0^)|*(l  -  \x\^)  (13) 

and 


WC2  =  2  * 


(14; 


The  first  term  WCl  takes  the  starting  value  for  the  water  content  computation  -  the  corresponding 
T3  value  -  and  reduces  it  exponentially  outward  from  the  z-axis  (i.e.,  from  x  =  0).  The  second 
term  WC2  provides  large  computational  values  near  the  z-axis  at  mid-level  heights  (i.e.,  near  H2) 
but  small  values  at  large  absolute  x  (i.e.,  near  the  periphery);  for  a  given  magnitude  of  x,  the  last 
factor  of  WC2  provides  a  skewness  in  x,  with  larger  (smaller)  numerical  values  at  positive 
(negative)  x. 

For  those  values  of  x  >  0.4  or  x  <  -0.6  units,  WC(x,0,z)  is  redefined  as  follows: 

WC(x,0;i)=Q  for  z<HB  (15) 

and 
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WCixM  =  WC ix,0^)  ^ 


forz>HB 


(16) 


Equation  15  is  important  since  it  both  defines  the  cloud  base  in  non-precipitating  areas  (as  HB) 
and  sets  the  horizontal  extent  of  surface  precipitation  (from  x  =  0.4  to  x  =  -0.6  units).  The  last 
two  factors  of  eq.  16  are  used  to  further  skew  WC(x,0,z)  in  a  positive  x  direction  at  middle  and 
high  levels. 

For  values  of  z  >  H6  (i.e.,  z  >  H/6),  the  full  3-D  representation  for  total  water  content  is 
given  as 

WC{x^;l)  =  WCixfi^)  *  YFAC  (17) 

Below  H6,  WC(x,0,z)  is  multiplied  by  an  additional  factor,  such  that 

WC(x^^)  =  WC(X,0;0  *  YFAC  *  (18) 

The  last  factor  of  eq.  18  provides  for  a  more  rapid  exponential  decrease  in  WC(x,o,z)  away  from 
the  x-axis  (in  the  y-direction),  and  results  in  more  circular  isopleths  of  WC  on  the  (x,y)  plane  at 
all  z  below  H6. 

Horizontal  Wind  WD.WS 

Generally,  the  spatial  structure  of  the  horizontal  wind  field  within  a  thunderstorm  is  highly 
complex,  varying  widely  from  storm  to  storm.  In  spite  of  this,  there  are  certain  wind  features 
that  relate  to  most  thunderstorm  occurrences  which  can  be  easily  simulated.  These  include  the 
strong,  divergent  surface  wind  outflow  below  the  storm's  downdraft  and,  at  mid-levels,  an  inner 
storm  structure  which  differs  considerably  from  that  some  distance  away.  For  the  latter,  the 
TSTORM  parameterization  assumes  that  the  centra)  core  of  the  storm  cell  (i.e.,  the  region  of 
strongest  up-  and  downdrafts)  offers  a  formidable  resistance  to  the  horizontal  wind.  Specifically, 
incoming  flow  is  both  diminished  and  deflected  on  its  approach  toward  the  inner  core;  once  past 
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this  core,  winds  converge  and  increase  before  exiting  the  storm  environment.  At  high  levels  near 
the  storm's  top,  well  away  from  the  strong  mid-level  core,  in-cloud  winds  are  relatively 
unchanged  from  the  outside  (environmental)  wind. 

Given  the  wind  speed  at  cloud  top  height  H  (VMAX).  TSTORM  defines  the 
environmental  wind  at  any  level  z  as 

(-1-68.:^)  (19) 

VE{z)  =  VMAX  *  e'  ^  ' 

For  a  grid  location  (x,y,z),  the  in-cloud  horizontal  wind  speed  (VC)  is  initially  computed  as 

KC  -  KE  *  RAD  (20) 

where  RAD,  the  horizontal  distance  from  the  storm  center  (x^^,  ,0)  to  the  location  (x.y)  at  level  z, 
is  defined 


RAD  =  (21) 

In  order  to  simulate  vertical  tilt  of  the  storm  with  respect  to  the  x-axis,  the  storm's  center 
x-coordinate  (x^.  )  is  defined  as 


-.16 

x^  = - z 

"  H6 


for  z  <  H6y 


(22) 


for  z>H6  and  z<  H-H6,  (23) 


and 


;C  =  I  2^8  ♦ 


+  0.08  for  Z  >  H-H6 


(24) 


where  H3  =  H/3.  Both  equations  (22)  and  (24)  provide  negative  vertical  tilts  wdth  respect  to  x, 
while  eq.  23  provides  a  positive  vertical  slope  from  z  =  H6  to  z  =  H  -  H6. 
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Given  values  of  VE  and  VC  at  the  grid  location  (x,y,z),  the  in-cloud  wind  speed  VC  is 
redefined  to  be  the  smaller  of  these  two  values.  Here,  VC  would  he  reset  to  VE  for  all  values 
of  RAD  greater  than  one.  The  wind  speed  VC  is  next  recomputed  in  terms  of  vertical  distance 
from  the  storm's  center  (H2),  that  is, 

KC  =  KC  +  {(VE  -VC)*-^*\z-  H2\]  (25) 

The  wind  speed  value  given  by  (25)  is  the  final  computational  value  for  all  z  >  H6.  For  levels 
z  <  H6,  the  strong  horizontal  wind  speeds  associated  with  the  cell's  downdraft  are  parameterized 
as 

VG  =  {VC  *  (26) 

For  a  grid  location  (x,y)  at  or  below  z  =  H6,  the  final  in-cloud  wind  speed  VC  is  defined  as  the 
larger  of  V'C  (from  eq.  25)  and  VG. 

Given  values  of  the  surface  wind  direction  (WDSFC)  and  the  surface  to  top  wind  veer 
(VEER),  the  environmental  wind  direction  (WDE)  at  any  level  z  is  defined  by 

WDE(z)  =  WDSFC  +  ^VEER  * 

In  those  cases  when  WDE  is  computed  to  be  greater  than  360  deg  (less  than  0  deg),  a  numerical 
value  of  360  is  subtracted  (added  to)  the  initial  value. 

At  all  levels  within  the  thunderstorm  cell,  the  wind  direction  (WDC)  is  set  equal  to  WDE 
provided  that  the  horizontal  distance  RAD  from  the  location  (x.y)  to  the  storm's  center  (x^,  ,0), 
as  computed  by  eq.  21,  is  greater  than  or  equal  to  one.  At  all  locations  where  RAD  is  less  than 
one,  TSTORM  parameterizes  the  wind  direction  to  differ  from  WDE  as  an  inverse  function  of 
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the  distance  to  the  storm's  inner  core.  This  dift’crencc  between  the  backiiround  wind  direction 


and  the  in-cloud  direction  WDC  at  a  location  (x,y.z)  is  designated  the  dedection  angle  ANG. 

To  determine  the  deflection  angle,  the  direction  vector  WDIV  from  the  storm  cell's  center 
(X(,  ,0)  to  the  location  (x,y)  is  first  specified  as  follows: 


WDIV  =  90. 

-  arctanly  /  x-x^) 

fur  x-x^ 

<  0 

WDIV  =  270. 

-  urctun(y  /  x-x^) 

for  x-x^ 

>  0 

WDIV  =  360. 

for  x-x^ 

=  0. 

y< 

0 

WDIV  =  999. 

for  x-Xf 

=  0, 

y  = 

0 

WDIV  =  180. 

for  x-x^ 

=  0, 

y> 

0 

(28) 


For  a  computed  value  of  WDIV  =  999,  the  corresponding  in-cloud  wind  direction  WDC  is  set 
to  999.  to  indicate  either  no  wind  (if  WS  =  0)  or  a  variable  wind  direction  (if  WS  0). 
Specifying  the  difference  between  the  background  wind  (WDE)  and  the  directiim  vector  WDIV 
as  WDIF,  the  magnitude  of  the  in-cloud  wind  deflection  (in  degrees)  is  initially  given  as 


lANGl  =  {(lWDIFl/90.)  -  INT (\ WDJF 1/90.)}  *  90.  (29) 

where  "INT"  signifies  the  integer  truncation  of  the  specified  value.  If  |ANG|  is  greater  than  45 
deg.  the  magnitude  of  the  deflection  angle  is  redefined  as  lANG]  =  90.  -  lANG|.  To  simulate 
vertical  variation  in  wind  deflection,  the  magnitude  |A.\’G|  is  moditled  as 


|ANG|  =  |AiVG|  ♦  ^  * 


(30) 


The  sign  of  the  deflection  angle  (positive  if  to  the  right  of  the  WDE  direction  vector,  negative 
if  toward  the  left)  is  computed  in  terms  of  the  polar  coordinates  (x'  .y'),  where 
x'  =  RAD  •  cos  (WDIF) 

and  (31) 

y’  =  RAD  *  sin  (WDIF) 
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For  x'  and  y'  both  greater  than  zero,  or  both  less  than  zero,  the  sign  of  the  dctlection  angle  is 
positive;  otherwise,  the  sign  of  ANG  is  negative. 

Given  the  deflection  angle  ANG.  the  in-cloud  wind  direction  for  a  location  (x,y)  at  a 
horizontal  distance  RAD  less  than  one  unit  from  the  storm's  center  is  simply 

VVDC  =  WDE  +  ANG  (32) 

Analogous  to  equation  27,  computed  values  of  WDC  greater  than  360  deg  (less  than  0  deg)  are 
brought  within  permissible  range  by  subtracting  (adding)  a  numerical  value  of  360  to  the  initial 
value. 

For  the  level  z  =  H6  and  below,  the  horizontal  wind  field  is  parameterized  to  diverge 
outward  in  all  directions  from  the  storm  center  (x^  ,0),  provided  that  the  horizontal  distance  from 
(x,y)  to  (Xj.  ,0)  is  less  than  or  equal  to 

RADDIV(z)  =  0.5  *  e*'  (33) 

Here,  wind  directions  at  any  location  (x,y)  are  specified  according  to  eq.  28  (i.e.,  WDC  = 
WDIV).  For  horizontal  distances  greater  than  RADDIV.  but  less  than  one  unit,  the  in-cloud  wind 
direction  WDC(x,y,z)  at  all  z  <  H6  remains  as  previously  computed  by  eq.  32. 

Dropsize  Distribution  DSD 

The  raindrop  size  distribution  at  a  user-specified  location  within  the  thunderstorm  cell  is 
computed  according  to  the  negative  exponential  relation  (Marshall  and  Palmer,  1948) 

DSDiNLOC.D)  =  (34) 

where  D  is  the  drop  diameter,  DSD(NLOC,D)  represents  the  number  of  drops  per  millimeter 
diameter  interval  and  per  cubic  meter  of  air,  is  a  constant  and  RR  is  the  site  rainfall  rate  in 
mm/hr.  The  constant  N^,  ,  which  corresponds  to  the  value  of  DSD(NLOC,D)  for  D  =  0,  is  set 
to  8000  mm'*  m'^  .  The  dropsize  distribution  is  calculated  for  21  diameters  D.  from  0  to  5  mm. 
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in  increments  of  6D  =  0.25  mm.  The  rainfall  rate  RR  is  determined  using  the  empirical  relation  that 
the  liquid  water  content  of  1  gm  ^  corresponds  approximately  to  the  precipitation  intensity  of  20 
mm/hr  (Blanchard,  1953). 

3.  TEST  CASE 

In  this  section,  vertical  and  horizontal  cross  sections,  and  raindrop  distributions, 

corresponding  to  the  TSTORM  simulation  listed  in  the  Appendix,  are  presented.  The  physical 

initialization  parameters  used  for  this  sample  case  are  as  follows: 

H  =  12  km  VMAX  =  15  m/s 

WMAX  =  9  m/s  WDSFC  =  90  deg 

VEER  =  +180  deg 

Based  on  these  numerical  values,  this  simulation  corresponds  to  a  thunderstorm  of  a  heavy 
intensity.  As  for  ail  TSTORM  simulations,  test  case  results  either  represent  the  structure  within 
an  isolated,  single-celled  thunderstorm  or,  the  structure  of  one  of  the  several  more-or-less 
independent  cells  found  within  the  common  multi-celled  thunderstorm. 

In  order  to  interpret  test  case  results  in  terms  of  spatial  and  temporal  dimensions,  the 
length  of  the  TSTORM  grid's  x-axis  (from  -1  to  +1  units)  is  assumed  to  correspond  to  a  distance 
of  20  km  and  a  time  of  one  hour.  These  values  are  chosen  to  agree  with  observations  of  Byers 
and  Braham  (1949),  who  determined  that  a  thunderstorm  cell  of  depth  H  =  12  km  has  a 
horizontal  extent  at  maturity  of  about  11  to  12  km,  and  that  the  average  time  duration  of  rainfall 
from  a  single  cell  of  a  storm  is  about  a  half  hour.  For  a  3-D  spatial  interpretation  of  test  results 
(i.e.,  X  corresponds  to  distance),  the  interval  of  (-0.7  <  x  <  0.5)  would  approximately  define  the 
spatial  structure  of  the  mature  thunderstorm  cell,  with  results  at  x  >  0.5  and  x  <  -0.7  units 
representing  cloudy,  non-precipitating  areas  on  the  periphery  of  the  storm.  Remembering  that  the 
storm's  translational  motion  (in  this  case,  20  km/hr)  is  assumed  by  TSTORM  to  be  in  the  positive 
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x-direction,  the  interval  from  x  =  0.4  to  x  =  -0.6  units  (which  defines  the  horizontal  extent  of 
surface  precipitation)  corresponds  to  exactly  a  half  hour. 

3.1  Vertical  Cross  Sections 

At  maturity,  the  vertical  velocity  structure  of  a  thunderstorm  cell  is  dominated  by  two 
separate  entities  -  the  updraft  and  downdraft.  The  updraft  results  from  the  strong  entrainment  of 
air  at  low  levels  into  the  thunderstorm  on  its  forward  flank,  while  the  downdraft  is  driven  by 
cold,  sinking  air  associated  with  the  descent  of  precipitation.  Measurements  reported  by  Byers 
and  Braham  (1949)  indicate  updraft  speeds  mostly  in  the  range  5  -  10  m/s  (these  averaged  over 
the  width  of  the  drafts)  with  maximum  values  most  common  at  middle  to  upper  levels  and, 
updraft  widths  typically  1  to  1.5  km  across.  The  vertical  cross  section  of  W  for  the  test 
simulation  (Figure  2)  indicates  good  agreement  with  these  observations,  showing  the  core  of  the 
updraft  to  be  at  (x  z)  coordinates  (0.3,  6.5)  with  a  magnitude  in  excess  of  6  m/s  over  a  1  to  1.5 
km  wide  horizontal  distance.  The  core  of  the  simulated  downdraft  is  found  at  about  5  to  5.5  km, 
with  a  magnitude  of  4.7  m/s;  these  results  are  in  agreement  with  Byers  and  Braham,  who  indicate 
that  downdrafts  are  less  intense  and  lower  than  updrafts.  In  Figure  2,  the  region  between  the 
updraft  and  downdraft  cores  is  one  of  very  strong  horizontal  shear  (i.e.,  turbulence).  Although 
constrained  by  the  TSTORM  parameterization  of  W  =  0  at  top  and  bottom,  careful  inspection  of 
Figure  2  reveals  that  the  updraft  tends  to  "shoot"  toward  the  middle  top  of  the  storm,  while  the 
downdraft  core  descends  conspicuously  toward  the  surface.  At  low  levels  (near  surface  to  2  km) 
below  the  updraft,  the  isotach  W  =  0,  which  represents  the  forward  outflow  boundary  of  the 
downdraft,  is  observed  to  have  a  pronounced  negative  vertical  tilt.  Finally,  at  values  of  absolute 
X  near  one,  vertical  velocities  depicted  in  Figure  2  are  quite  small;  such  benign  values  represent 
pre-  and  post-storm  conditions  or,  vertical  motion  sufficiently  removed  from  the  storm's 
periphery. 
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The  temperature  anomaly  vertical  cross  section  (Figure  3)  is  dominated  by  two  opposing 
cores.  Within  the  storm's  forw/ard  sector,  a  positive  anomaly  core  of  about  2.5  C  is  located  a  few 
kilometers  below'  the  updraft  maximum.  Related  to  the  storm's  downdraft,  cold  temperature 
anomalies  of  magnitude  1  to  2  C  are  found  at  mid-levels,  increasing  appreciably  tow  ard  the 
surface,  with  a  maximum  negative  anomaly  in  excess  of  7  C  at  the  surface  below  the  downdraft. 
The  0  C  isotherm  which  underlies  the  updraft  in  a  wedge-like  fashion  corresponds  to  the  leading 
edge  of  the  cold  air  outflow  beneath  the  storm.  Note  that  the  temperature  anomaly  gradient  at 
the  surface  is  very  strong  from  the  time  of  the  arrival  of  the  outflow  (x  =  0.4)  to  the  time  of  the 
storm's  apex  (x  =  0).  Once  beyond  the  storm’s  center,  the  surface  temperature  anomaly  rises 
steadily  before  leveling  off  at  about  -1.5  C  for  x  <  -0.5  .  These  negative  surface  temperature 
anomalies  at  large  negative  x  simulate  the  cool  air  left  in  the  storm's  wake.  On  the  forward 
perimeter  of  the  storm,  at  the  pre-storm  stage  (i.e.,  near  x  =  1),  TSTORM  produces  little  if  any 
anomaly  from  the  outside  environment  at  any  level. 

The  vertical  cross  section  of  the  pressure  perturbation  (Figure  4)  shows  two  opposing 
anomaly  cores  at  mid-level;  a  weaker  negative  Pp  core  of  about  -0.7  mb  within  the  downdraft  and 
a  stronger  positive  Pp  core  exceeding  1.2  mb  within  the  updraft.  At  the  surface  and  lower  levels, 
the  TSTORM  test  case  exhibits  a  very  slight  pressure  reduction  below  the  updraft  and  a 
prominent  pressure  dome  of  about  +0.7  mb  below  the  downdraft  at  the  center  of  the  storm. 
Inspection  of  tabular  results  of  Figure  4  reveals  that  the  surface  pressure  dome  anomaly  initially 
decreases  after  the  storm's  apex,  then  levels  off  to  about  +0.15  mb  at  large  negative  x  (i.e.,  at  x 
<  -  0.5). 

As  given  by  the  vertical  cross  section  of  total  water  content  (Figure  5),  the  shape  of  the 
test  simulation  thunderstorm  cell  is  basically  cylindrical.  However,  on  close  inspection,  one  can 
discern  a  skewness  in  total  water  content  at  middle  to  upper  levels  in  the  direction  of  storm 
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Figure  4.  Test  case  vertical  cross  section  (at  y  =  0)  of  pressure  perturbation  P 
Plot  contours  range  from  -0.6  to  1.2  mb,  at  intervals  of  0.2  mb. 
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Figure  5.  Test  case  vertical  cross  section  (at  y  =  0)  of  total  water  content  WC. 
Plot  contours  range  from  0.5  to  7.0  g  m'^ ,  at  intervals  of  0.5  g  m'^  . 


motion  (i.e.,  in  the  +x  direction),  with  significantly  higher  values  of  WC  ahead  of  the  storm  (x 
>  0.8)  than  in  the  rear  (x  <  -0.8).  Within  the  thunderstorm  cell,  a  large  concentration  of  water 
content  (>  5.5  g  m'^  )  is  found  at  middle  levels,  covering  the  inner  storm  region  between  the 
updraft  and  downdraft  cores.  Below  3  km,  the  storm's  highest  total  water  content  becomes 
concentrated  within  a  narrow  band,  forming  a  rain  shaft  which  drops  to  the  surface  at  the  storm's 
center  (x  =  0).  Based  on  Blanchard's  empirical  relation,  the  maximum  value  WC  =  7  gm’^  at  the 
storm's  apex  corresponds  to  a  (3  minute  average)  rainfall  rate  of  140  mm/hr.'  This  simulated 
thunderstorm  would  produce  a  total  of  about  21  mm  (0.83  in)  of  rain  in  about  half  an  hour  (i.e., 
during  the  time  -0.6  <  x  <  0.4)  at  a  site  directly  in  the  storm's  path. 

In  the  test  simulation,  the  vertical  wind  structures  during  pre-  and  post-storm  conditions 
are  nearly  identical.  At  these  times,  horizontal  winds  are  observed  to  veer  steadily  with  height 
(from  surface  easterly  to  upper  level  westerlies),  with  speed  increasing  gradually  from  the  surface 
to  mid-levels,  then  more  rapidly  to  the  storm's  top  (Figures  6  and  7).  Within  the  thunderstorm 
cell  close  to  its  vertical  axis  (i.e.,  near  x  =  0),  the  wind  structure  is  much  more  varied.  At  the 
top,  westerly  winds  of  15  m/s  decrease  rapidly  toward  a  mid-level  core  of  very  light  winds,  with 
an  absolute  minimum  of  0  m/s  at  (x,z)  =  (0,6).  Below  this  inner  storm  minimum,  winds  first 
increase  slowly  toward  the  surface,  then  very  rapidly  once  inside  the  surface  cold  air  downdraft. 
Within  the  confines  of  the  cold  air  outflow,  wind  directions  are  divergent  from  the  storm's  center. 
The  arrival  of  the  surface  "gust  front"  is  marked  by  a  sudden  shift  in  wind  direction  (from  90  deg 
at  X  =  0.6  to  270  deg  at  x  =  0.5)  and  a  rapid  increase  in  wind  speed  from  3  to  10  m/s  in  about 
10  minutes.  Once  the  apex  of  the  storm  passes,  surface  horizontal  wind  speeds  decrease  rapidly, 
falling  below  3  m/s  at  x  =  -0.4. 

‘For  this  and  any  simulation,  the  storm  maximum  rainfall  rate  can  be  found  by  multiplying 

the  storm  depth  H  (km)  by  11.7  mm  hr  '  km  ‘. 
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Figure  7.  Test  case  vertical  cross  section  (at  y  =  0)  of  horizontal  wind.  Plot  scale  represents  a 
westerly  wind  of  10  m/s.  Table  gives  wind  direction  WD  at  actual  grid  point  resolution. 
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3.2  Horizontal  Cross  Sections 

For  the  test  case  simulation,  plots  of  horizontal  (x,y)  cross  sections  were  made  for  all 
computed  meteorological  fields  for  both  the  surface  (except  z  =  0.5  km  for  W)  and  the  mid-storm 
level  (z  =  6  km).  In  order  to  easily  interpret  horizontal  dimensions  of  thunderstorm  features,  the 
length  of  the  TSTORM  grid  y-axis  is  assumed  equal  to  that  of  the  x-axis,  that  is,  20  km. 

Given  the  bottom  boundary  condition  W  =  0,  values  of  the  vertical  velocity  just  above  the 
surface  (at  z  =  0.5  km)  are  predictably  small  (Figure  8a);  nonetheless,  the  thunderstorm  cell's 
almost  circular  downdraft  is  clearly  evident  near  the  center  of  the  TSTORM  grid.  At  mid-level 
(Figure  8b),  a  iwo-lobed,  updraft-downdraft  pattern  is  well  defined.  Here,  both  the  updraft  and 
downdraft  exhibit  elliptical  structure,  with  the  outer  contours  of  each  appreciably  elongated  in  the 
y-direction. 

The  horizontal  cross  sections  of  temperature  anomaly  at  the  surface  and  z  =  6  km  are 
quite  distinct  (Figure  9).  At  mid-level,  two  opposing  cores  of  T^  of  similar  magnitude  (about  2 
C)  are  evident;  much  like  the  vertical  velocity,  the  shape  of  the  positive  and  negative  T,  patterns 
are  elliptical,  being  elongated  in  the  y-direction.  The  (x,y)  cross  section  plot  at  the  surface  is 
dominated  by  the  very  prominent  temperature  deficit  a.ssociated  with  the  cold  air  outflow  beneath 
the  downdraft.  Contours  of  this  feature  are  observed  to  be  slightly  elliptical  within  0.4  units  of 
the  center.  The  open  contours  at  x  <  -0.6  represent  cooler  surface  conditions  in  the  wake  of  the 
storm,  that  is,  post-storm  locations  where  the  temperature  has  not  recovered  to  pre-storm  levels. 

In  general,  the  horizontal  cross  sections  of  the  pressure  perturbation  appear  quite  similar 
to  those  of  the  vertical  velocity.  At  the  surface  (Figure  10a),  the  closed  contours  near  the  grid 
center  describe  the  spatial  structure  of  the  storm's  pressure  dome  beneath  the  downdraft.  At 
mid-level  (Figure  10b),  a  two-lobed,  elliptical  contour  structure  is  evident,  with  the  positive  Pp 
core  (related  to  the  updraft)  about  double  the  magnitude  of  the  weaker  negative  Pp  core. 
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Figure  8.  Test  case  horizontal  cross  sections  of  vertical  velocity  W.  a)  at  z  =  0.5  km,  b)  at 
z  ■=  6  km.  Plot  contours  range  from  *1.0  to  0.0  m/s  in  a)  and  from  -4.5  to  6.0  m/s  in  b), 
both  in  intervals  of  0.5  m^s. 
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Figure  9.  Test  case  horizontal  cross  sections  of  temperature  anomaly  T,.  a)  at  z  =  0  km,  b)  at 
z  =  6  km.  Plot  contours  range  from  -7.0  to  0.0  C  in  a)  and  from  -2.0  to  2.0  C  in  b),  both 
in  intervals  of  0.5  C. 
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The  surface  horizontal  cross  section  of  the  total  water  content  (Figure  11a)  delineates  the 
spatial  structure  of  the  precipitation  field.  Here,  the  smallest  contour  (0.5  g  m  ^  )  corresponds 
to  a  moderate  rainfall  rate  of  10  mm/hr;  the  area  within  this  contour  is  quite  small,  being  about 
40  km^.  As  to  be  expected,  the  very  intense  rainfall  rate  (RR  >  100  mm/hr)  of  the  simulated 
thunderstorm  is  concentrated  over  a  very  small  area  (several  km  )  at  the  storm's  center. 

The  horizontal  structure  of  the  thunderstorm  cloud  at  mid-level  is  given  in  Figure  11b. 
Here,  the  shape  of  the  WC  contours  are  elliptical,  with  the  x-axis  being  the  major  axis.  The  inner 
contours  of  high  total  water  content  are  markedly  flattened;  note  the  lateral  extension  in  the 
negative  x-direction  from  the  core  maximum  at  x  =  0.2  units.  At  this  level  (z  =  6  km),  the  con¬ 
tours  of  WC  are  noticeably  skewed  toward  the  right  (i.e.,  in  the  +x  direction).  Given  the 
assumed  (x,y)  spatial  dimensions,  the  thunderstorm  cloud  at  this  level  extends  over  a  horizontal 
area  of  a  few  hundred  square  kilometers. 

As  formulated  in  TSTORM,  contours  of  horizontal  wind  speed  are  essentially  circular. 
In  Figure  12a,  the  outer  contour  (WS  =  3  m/s)  represents  the  outward  extent  of  winds  (i.e.,  the 
gust  front)  beneath  the  storm;  in  terms  of  assumed  horizontal  dimensions  for  this  simulation,  this 
boundary  occurs  at  about  4  km  from  the  storm's  center.  At  mid-level  (Figure  12b),  horizontal 
winds  on  the  storm's  perimeter  (in  this  simulation,  of  speed  6.5  m/s)  decrease  steadily  toward  the 
center  core  of  the  storm. 

The  surface  wind  vectors  of  Figure  13a  indicate  that  within  the  confines  of  the  gust  front, 
winds  blow  radially  outward  from  the  storm  center,  with  increasing  speed  toward  the  center;  out¬ 
side  of  the  gust  front,  winds  are  uniform  (viz.,  easterly  at  2.8  m/s).  At  mid-storm  level  (Figure 
13b),  horizontal  winds  are  deflected  and  slow  down  as  they  approach  the  storm  center;  once  past 
the  inner  core,  winds  converge  and  speed  up,  eventually  exiting  the  storm  environment  with  a 
wind  speed  and  direction  equal  to  that  of  the  environmental  wind  (viz.,  217  deg  and  6.5  m/s). 
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Figure  12.  Test  case  horizontal  cross  sections  of  wind  speed  WS.  a)  at  z  =  0  km,  b)  at 
z  =  6  km.  Plot  contours  range  from  3.0  to  10.0  m/s  in  a)  and  from  0.5  to  6.0  m/s 
in  b),  both  in  intervals  of  0.5  m/s. 
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Figure  13.  Test  case  horizontal  cross  sections  of  wind  vectors,  a)  at  z  =  0  km,  b)  at  z  =  6  km. 
Plot  scale  represents  a  westerly  wind  of  10  m/s.  Note:  at  coordinates  (0,0,0),  the  wind 
direction  is  variable;  at  (0,0,6),  the  horizontal  wind  is  zero. 
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3.3  Dropsize  Distribution 

In  the  test  simulation,  Marshall-  Palmer  dropsize  distributions  were  computed  at  four  sites 
within  the  lower  thunderstorm  cloud  structure.  Figure  14  presents  these  dropsize  distribution 
results  in  both  graphical  and  tabular  form.  For  interpretation  of  results,  the  actual  number  of 
drops  per  cubic  meter,  of  diameter  between  D  and  D  +  6D  (where  6D  is  preset  to  0.25  mm),  is 
obtained  as  the  product  of  the  distribution  function  DSD  and  6D;  for  example,  at  site  a),  the 
number  of  drops  of  diameter  3.25  to  3.5  mm  in  a  cubic  meter  of  air  would  be  68  *  0.25,  or  17 
drops. 

At  any  particular  site,  and  for  a  specified  drop  diameter,  the  actual  dropsize  distribution 
is  determined  by  the  rainfall  rate  RR  which,  in  turn,  is  a  direct  function  of  the  total  water  content. 
Generally,  for  a  fixed  drop  diameter,  the  number  of  drops  per  cubic  meter  of  air  is  in  direct 
proportion  to  the  rainfall  rate  -  the  higher  the  rate,  the  more  drops.  Of  the  sites  chosen  for  this 
simulation,  the  dropsize  distribution  is  largest  at  site  a),  which  is  located  just  above  the  surface 
in  the  rain  shaft,  and  least  at  site  d),  located  slightly  above  the  ground  in  lightly  descending  air 
at  X  =  -0.6  units.  Careful  inspection  of  Figure  14  indicates  that,  while  the  actual  difference  in 
the  number  of  drops  between  these  two  sites  is  largest  at  small  diameters  (excepting  D  =  0),  the 
percent  difference  is  largest  at  large  diameters,  where  there  are  far  fewer  drops  per  cubic  meter. 
For  all  sites  a)  through  d),  the  distribution  function  DSD  (in  mm’ m'^ )  decreases  exponentially 
with  drop  diameter,  from  several  thousand  at  small  drop  diameters  to  about  10  or  less  at 
diameters  near  5  mm. 
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Figure  14.  Test  case  dropsize  distributions  at  selected  cloud  locations,  a)  at  (0,0, 0.5) 
b)  at  (-0.2, 0,2.0),  c)  at  (0.4,0,2.0)  and  d)  at  (-0.6,0,!. 0).  The  rainfall  rate  RR 
(mm/hr)  is  also  given  for  each  location. 
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APPENDIX 


TSTORM  COMPUTER  LISTING 


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


TSTORM 

A  FEATURE  MODEL  OF  A  THUNDERSTORM  CELL 
CCCCCCCCCCCCCCCCCCC  MAIN  PROGRAM  C CC C CC C C C C C C C C C C C C C C C CC C C C C C C C C C C C C C 


INPUTS  VARIABLE 
NLVL 
H 

MMAX 

VMAX 

WDSFC 

VEER 

NLOC 

XOSD(NLOC) 
YOSO<NLOC) 
ZDSD ( NLOC  ) 


DESCRIPTION 
NUMBER  OF  LEVELS 
THUNDERSTORM  HEIGHT  [E 

MAX.  VERT.  VELOCITY 
WIND  SPEED  AT  LEVEL  H 
SFC.  WIND  DIRECTION 
WIND  VEER,  SFC.  TO  H 
NO.  OF  SITES  FOR  DROPSIZE  DSTR. 

''•COORDINATE  FOR  NLOC  SITE  [-1 

Y-COORDINATE  FOR  NLOC  SITE  {-1 

Z-COORDINATE  FOR  NLOC  SITE  [0. 


]  KM 


(VALID  RANGE, UNITS) 
13-31 
.14.6,16 
16.  M/S 
30. 

360  . 

360. 

>«  1 
.-0.9. . 

,-0.9.. 

0 . 6  ...  4 


6.6 

6 

10  . 
0  .  ■ 


M/  S 
DEC 
DEG 


.0.9. 

.0.9. 

.6.6. 


] 


] 

] 

KM 


OUTPUTS  VARIABLE 

W(Z1  ,21  ,NLVL) 
TA(21  ,21  ,NLVL) 
PP(21  ,21  ,NLVL) 
WC( 21  .21  ,  NLVL  ) 
WS{21  .21  .NLVL) 
WD(21  ,21  .NLVL) 
00(21) 

DSD ( NLOC  .  21  ) 

RR ( NLOC  ) 


DESCRIPTION 

VERTICAL  VELOCITY  ARRAY 
TEMPERATURE  ANOMALY  ARRAY 
PRESSURE  perturbation  AHRAY 
TOTAL  WATER  CONTENT  ARRAY 
WIND  SPEED  ARRAY 
WIND  DIRECTION  ARRAY 
DROP  diameter  array  [0 

DROPSIZE  DISTRIBUTION  ARRAY 
RAINFALL  RATE  AT  NLOC  SITE 


(VALID  RANGE. UNITS) 
M/S 
C 

MB 

C/M3 

M/S 

DEG 

.  .  0 . 25 .  .  . 4 . 76 , 6 . ]  MM 
MM/HR 


THE  TSTORM  ALGORITHM  IS  A  FEATURE  MODEL  WHICH  PROVIDES  A  SIMPLE 
PARAMETRIC  R E P R E S F N T A T I  ON  OF  A  THUNDERSTORM  CELL.  IT  IS  DESIGNED 
FOR  USE  IN  SIMULATIONS  REQUIRING  AN  IDEALIZED  MAGNITUDE  AND  SPATIAL 
DISTRIBUTION  OF  METEOROLOGICAL  PARAMETERS  MITHIN  A  THUNDERSTORM. 
TSTORM  OUTPUT  MAY  BE  SUPERIMPOSED  ON  AN  INITIAL  DATA  FIELD  TO 
PROVIDE  ADDED  SUBGRID  INFORMATION  OR  TIME  VARIABILITY.  OR  MAY  BE 
USED  AS  A  UNIQUE  ENVIRONMENTAL  SOURCE  FOR  THE  HOST  SIMULATION. 


THE  FOLLOWING  PROGRAM.  WRITTEN  IN  ANSI  FORTRAN  77.  IS  A  DEMONSTRA¬ 
TION  DRIVER  FOR  THE  TSTORM  SUBROUTINE.  BASED  ON  THE  PHYSICAL  INPUT 
PARAMETERS,  THIS  SIMULATION  CORRESPONDS  TO  A  THUNDERSTORM  OF  A 
HEAVY  INTENSITY.  AFTER  RETURN  FROM  THE  TSTORM  SUBROUTINE,  THE 
DRIVER  PROGRAM  OUTPUTS  (IN  TABULAR  FORM)  SEVERAL  VERTICAL  CROSS 
SECTIONS  AND  DROPSIZE  DISTRIBUTIONS  FOR  THE  SIMULATED  THUNDERSTORM. 

START  DEMO  PROGRAM 

PROGRAM  TSMODL 

ENTER  ARRAY  PARAMETERS 

PARAMETER(NLVL«26) 

PARAMETER(NL0C>4) 

CHARACTER'10  TITLE1 (6) 

CMARACTER’*8  TITLE2(1) 

DIMENSION  W(21,21,NLVL),TA(21.21,NLVL),PP(21,21,NLVL) 

DIMENSION  WC(21,21,NLVL).WS(21.21.NLVL).WC(21.21.NLVL) 

DIMENSION  00(21 ).OSO(NLOC. 21 }.RR(NLOC) 

DIMENSION  XDSO(NLOC}.YDSO(NLOC).ZDSD(NLOC) 

DATA  TITLE1 (1 )/ '  W  (M/S)  ' / 

DATA  TITLE1 (2)/'  TA  (C)  '/ 
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DATA  TITLE1 (3)/-  PP  (MB)  '  ! 

DATA  T1TLE1(4)/-  «tC  (G/M3)'/ 

DATA  T I TLE1  (  5  )  /  '  NS  (M/ S  )  '  / 

DATA  TITLE1(6)/-  NO  (OEG)  '/ 

DATA  T1TLE2(1  )/  '  Y  DSPL«'  / 

ENTER  OSD  LOCATIONS 

DATA  XDSD/0 . 0 . 0 .4 . -0 . 2 . -0 . 6/ 

DATA  YDSD/0 . 0 , 0 . 0 . 0 . 0 , 0 . 0/ 

DATA  2DSD/0 . 5. 2. 0,2. 0,1.0/ 

YPLT  «  0.0 

ENTER  PHYSICAL  PARAMETERS 

H  «  1  2  . 

NMAX  -  9. 

VMAX  -  15. 

WOSFC  >  90. 

VEER  c  180. 

0PEN(7.FILE«'0UTFIL'  ,STATUS='NEN'  ) 

CALL  TSTORM  SUBROUTINE. 

CALL  TSTORM(NLVL.H.NMAX.VMAX.NDSFC.VEER.NLOC.XDSD.YDSD.ZOSD. 
•N,TA,PP,NC,WS,ND.OD.OSD.RR) 

NRITE  OUTPUT  FILES. 


C 

C 

C 

C 

c 


c 

c 

c 

c 

c 

c 


YR  =  (YPLT  -f  1.1)  •  10. 

J  e  INT ( YR) 

00  100  NP  «  1,6 

NRITE(7,200)  TITLE1(NP),TITLE2(1).YPLT 
200  FORMAT (A1 0  ,  A8  .  F4  .  1  ) 

DO  300  K  -  1 , NL VL 

KL  «  1  ♦  (NLVL-K) 

IF{NP.EQ.1)  NRITE{  7,400)  ( N ( I , J , K L )  ,  I  «  1  ,  2  1  ) 
IF{NP,EQ.2)  NRITE(  7,400)  ( T a ( I  ,  J  ,  K L  )  ,  I  ■ 1  ,  2 1  ) 

IF(NP.EQ.3)  NRITE(  7,400  )  ( P P (  I  ,  J  .  K L  )  ,  1 » 1  ,  2 1  ) 
IF(NP.EQ.4)  NRITE(  7.400)  { N C ( I , J  .  K L  )  ,  I  =  1  ,  2 1  ) 

IF(NP.EQ.5)  write (7, 400)  (NS(I.J,KL),1«1,21) 
IF(NP.EQ.6)  WRITE(  7,500)  ( WD (  I  ,  J  ,  K L  )  ,  I  » 1  ,  2 1  ) 

400  FORMAT ( 21 F5  .  2  ) 

500  FORMAT (21F5  .  1  ) 

300  CONTINUE 

100  CONTINUE 

DO  600  NL  -  1 , NLOC 

WRITE(7,700)  XOSO(NL),YOSO(NL),ZDSD(NL).RR(NL) 

700  FORMAT(-  X»-,F6.1,'  Y«',F5.1,-  2='.F5.1,'KM  RR=', 

•  F6 . 1 , ' MM/HR ' ) 

WRITE(7,600) 

800  F0RMAT('  0(MM)  OSD  {1/(M«*3  MM))') 

DO  900  ND  «  1,21 
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1000 

c 

900 

C 

600 

C 


WRI TE (  7  .  1  000  )  00  ( NO  )  .  OSO ( NL  .  ND  ) 
FORMAT (F5 . 2 , 6X  ,  F7  .  1  ) 

CONTINUE 

CONTINUE 

STOP 


ENO  OEMO  PROGRAM 


ENO 

CCCCCCCCCCCCCCCC  START  SUBROUTINE  TSTORM  C C CC C C C CC C CC CCC CC C C C C C C C CC C C 


SUBROUTINE  TSTORM(NLVL.H. MMAX . VMAX .WDSFC.VEER.NLOC.XOSO.YOSD.ZOSD. 
•N.TA.PP.WC.BS.WD.OO.OSD.RR) 


INPUTS  VARIABLE  DESCRIPTION 

NLVL  NUMBER  OF  LEVELS 

H  THUNDERSTORM  HEIGHT  [6 

MMAX  MAX.  VERT.  VELOCITY 

VMAX  WIND  SPEED  AT  LEVEL  H 

NOSFC  SFC.  WIND  DIRECTION 

VEER  MIND  VEER  ,  SFC .  TO  H 

NLOC  NO.  OF  SITES  FOR  DROPSIZE  DSTR. 

XOSO(NLOC)  X-COOROINATE  FOR  NLOC  SITE  ['1 

YOSD(NLOC)  Y-COOROINATE  FOR  NLOC  SITE  [-1 

ZOSO(NLOC)  Z-COOROINATE  FOR  NLOC  SITE  [0. 


(VALID  RANGE. UNITS) 
13-31 

..6.5.  ..14.6,16. 3  KM 
6  .  -  15.  M/S 
10.-  30 .  M/ S 
0.-  360.  DEG 
360.  DEG 
>»  1 

.-0.9. ..0.9.1.] 
.-0.9.  ..0.9,1. 3 
0.5. .  .4.6,6.  3  KM 


OUTPUTS  VARIABLE 

M{21  ,21  ,NLVL) 
TA{21  ,21  ,NLVL) 
PP (21, 21, NLVL) 
MC (21, 21, NLVL) 
MS( 21  ,  21  ,NLVL  ) 
M0( 21  ,  21  ,NLVL  ) 
00(21) 

OSD(NLOC,21  ) 
RR(NLOC) 


DESCRIPTION 

VERTICAL  VELOCITY  ARRAY 
TEMPERATURE  ANOMALY  ARRAY 
PRESSURE  PERTURBATION  AMAY 
TOTAL  MATER  CONTENT  ARRAY 
MIND  SPEED  ARRAY 
MIND  DIRECTION  ARRAY 
DROP  DIAMETER  ARRAY  [0 

DROPSIZE  DISTRIBUTION  ARRAY 
RAINFALL  RATE  AT  NLOC  SITE 


(VALID  RANGE, UNITS) 
M/S 
C 

MB 

G/M3 

M/S 

DEC 

. ,0.26. . .4.76,6. 3  MM 

MM/HR 


DIMENSION  M(21,21,NLVL).TA(21.21,NLVL).PP(21,21.NLVL) 
DIMENSION  MC(21,21,NLVL),MS(21,21,NLVL),MD(21,21,NLVL) 
DIMENSION  DO(21),OSO(NLOC,21),RR(NLOC) 

DIMENSION  XDSO(NLOC),YDSO(NLOC),ZDSD(NLOC) 


P I  -  4 . "ATAN ( 1  .  ) 
OEGRAO  -  180. /PI 


HB  B  H/ 1 2  . 

H6  B  H/6. 

H3  B  h/3. 

H2  B  H/2. 

HT  B  H  -  H/1 2  . 

DO  100  K  B  1 , NLVL 

Z  B  -0.5  ♦  0.5*K 

00  200  I  B  1,21 

RX  B  -1.1  ♦  0.1*1 


CALCULATE  THE  VERTICAL  VELOCITY. 

RXOPL  B  RX  -  (  0.432  •  ( ( H2 - Z  ) / H 2  )  • • 3  ) 

Ml  B  EXP( -PI*RX*RX) 


ooo  ooo  o  o  noo  ooo  o  ooo 


c 

c 


H2  >  SIN(PI*RXOPL) 

N3  >  «1*N2 

«4  •  (4-lfMAX/M)«(Z  -  2*2/H) 
NS  ■  (1  Z*RX/H) 

NCAL  >  N3*H4-HS 


CALCULATE  THE  TEMPERATURE  ANOMALY. 

T1  ■  EXP{-PI*RX*RX)  •  (0.3/H2)  •  (M-Z)**2 
T2  ■  1  ♦  { ( Z* Z*ABS(RX  )  ) / (M*H/24 .  )  ) 

T3  ■  EXP(RX)  •  RX  •  (  4.32  *  ( ( H2 - 2 ) / H 2 ) • • 3  ) 

TACAL  ■  N3«T1«T2 

IF(RX.LE.O.  .AND.  2.LE.H2)  TACAL  *  TACAL  ♦  T3 

CALCULATE  THE  PRESSURE  PERTURBATION. 

1F(2.LE.HB  .OR.  Z.GE.HT)  PPCAL  «  0.1*  (NCAL-TACAL) 

I F ( 2 . GT . HB  .  AND  .  2  .  LT  .  HT  ) 

PPCAL  »  0.1*  (NCAL-TACAL)  •  (2  -2*Z/H) 

CALCULATE  THE  TOTAL  NATER  CONTENT  (LIQUID  «  SOLID) 

NCI  ■  ABS(TACAL)  *  (  1  -  { AB S ( R X  )  * • 0 . 3 3 3  ) ) 

NC2  «  2  *  (2  -2*Z/H)  *  EXP  (  -PI  *((RX-  0 . 1 2 • 2 / H 2  )  • • 2 ) ) 
NCCAL  >  NCI  *  NC2 

IF ( (RX . GE . 0 . 4  .OR.  R X . L T . - 0 . 6  )  .  A N D  .  2.LE.HB)  NCCAL  -  0.0 
I F ( ( RX . GE . 0 . 4  .OR.  RX  .  L T  .  •  0 . 6  )  .  A N 0  .  2.GT.HB) 

NCCAL  «  NCCAL  •  EXP(-RX*RX)  *  (1  ♦  2»RX/M) 

00  300  J  -  1,21 

RY  «  -1.1  ♦  0.1‘J 
RYFAC  *  (1.0  -  ABS(RY)  )»*2 

SPECIFY  NCAL , TACAL . PPCAL  AND  NCCAL  AT  (I.J.K)  COORDINATES. 

N( I  . J  ,  K  )  «  NCAL  •  RYFAC 
TA( I  ,  J  ,K)  -  TACAL  •  RYFAC 

PP(I.J,K)  »  PPCAL  •  RYFAC 

NC(I,J,K)  ■  NCCAL  •  RYFAC 

IF(Z.LT.H6)  NC(I.J,K)  «  NC(I.J.K)  *  ( R Y F A C • • ( H6  -  2  )  ) 

CALCULATE  THE  HORIZONTAL  NINO  SPEED  AND  DIRECTION. 

VE  «  VMAX  •  EXP(-1.6e  •  (H-2)/H  ) 

C 

NOE  -  NDSFC  ♦  (VEER  *  S  I  N ( ( P  I  • 2  )  /  (  2  •  M  )  )  ) 

IF (NOE  .  GT  .  360  .  )  NDE  «  NOE  -  360. 

IF  (NDE  .  LE  .  0  .  )  NOE  ■  NOE  360. 

C 

IF{Z.LE.H6)  RXH  ■  RX  ♦  (  (0.16/M6)  *  2  ) 

IF(Z.GT.H6  .AND.  2.LT.(H-H6)) 

•  RXH  ■  RX  ♦  (  (0.16/H3)  •  (H6-Z)  )  ♦  0.16 

IF (2 . GE . (H-H6 )  )  RXH  ■  RX  ♦  (  (0.08/H6)  •  (Z-H)  )  -  0.08 
C 

RAO  •  SQRT(  (RXH«RXH)->(RY«RY)  ) 

C 

VC  ■  VE  •  RAO 
VC  >  AMIN1 (VE.VC) 

VC  ■  VC  ■»  (  ( VE  -  VC  )  •  ( 1 /H2)»ABS(  2 -H2  )  ) 

C 

IF(RA0.GT.  0.999}  GO  TO  400 
C 

C  CALCULATE  A  DEFLECTION  ANGLE  NITHIN  THE  STORM. 
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c 

I F ( RXH . LT  .  0  .  )  WDJV  >  90.  -  (DEGRAD  *  A T A N ( R Y / R X H ) ) 

IF ( RXH  .  GT  . 0 .  )  WDIV  .  270.  -  (DEGRAD  •  A T A h ( R Y / R X H  )  ) 
IF(RXM.EQ.O.  .AND.  RY.LT.O.)  NDIV  >  360. 
IF(RXH.EQ.O.  .AND.  RY.EQ.O.)  NO  1 V  .  999. 
IF(RXH.EQ.O.  .AND.  RY.CT.O.)  NDIV  «  160. 

C 

1F(ND1V.EQ.999  )  GO  TO  SOO 
C 

NOIF  -  NOE  -  NDIV 

C 

XROT  >  RAO  ■  COS ( NO  1 F /DEGRAD  ) 

YROT  -  RAO  •  SIN(NDIF/0EGRAD  ) 


SIGN  >  0  . 


I F (XROT . GT . 0 . 

.  AND  . 

YROT .GT  . 0  .  ) 

SIGN  E 

♦  1 

IF(XROT.LT.O. 

.  AND  . 

YROT .GT . 0  .  ) 

SIGN  E 

-  1 

IF(XROT .LT .0. 

.  AND  . 

YROT . LT . 0 . ) 

SIGN  E 

♦  1 

IF(XROT .GT .0. 

.  AND  . 

YROT . LT . 0 . ) 

SIGN  E 

-  1 

BETA  -  ABS(ND1F ) /90  . 

GAMMA  -  INT(BETA) 

ANC  >  ( beta  -  GAMMA ) *90 . 

I  F ( ANG  .  GT  . 45 .  ]  ANG  ^  90.  -  ANG 

ANG  s  ANG  •  SIGN 
C 

ANG  «  ANG  •  EXP (  -RAO'RAD  )  *  {H3/H)*(Z-  Z*Z/H) 

C 

NOC  -  NOE  «  ANG 

IF ( NDC  .  GT . 360 .  )  NOC  «  NOC  -  360. 

IF (NOC  .  LE  . 0 .  )  NOC  «  NDC  «  360. 

GO  TO  600 
C 

400  NOC  >  NOE 

GO  TO  600 
C 

SOO  NOC  «  999. 

C 

600  : F ( 2  .  GT . H6 )  GO  TO  700 

C 

VG  ■  (  VC  ♦  EXP(H6-Z))  •  EXP(-PI«RAD) 

c 

VC  ■  AMAXI  ( VC  ,  VG  ) 

C 

RADO I V  «  0.5  •  EXP(-2) 

C 

IF (RAD . LE . RAOOI V  )  NDC  «  NDIV 
C 

700  NS(I.J,K)  «  VC 
NO ( I  , J , K  )  -  NDC 
C 

300  CONTINUE 
C 

200  CONTINUE 
C 

100  CONTINUE 

CALCULATE  THE  DROPSIZE  DISTRIBUTION  AT  SELECTED  LOCATIONS. 
DO  BOO  NL  -  1 , NLOC 


XR 

-  (  XOSD(NL) 

♦  1.1) 

•  10 

I  ■ 

INT(XR) 

YR 

■  (  YDSO(NL) 

♦  1.1) 

•  10 

J  ■ 

INT{YR) 

ZR 

>  (  ZOSO(NL) 

♦  0.5) 

•  2  . 

K  E 

INT(ZR) 
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c 

c 

RR(NL)  ■  20.0  *  «C(I.J.K) 

ALPHA  «  4  .1  •  ( 1 /RR(NL  )**0 . 21  ) 

C 

00  900  NO  -  1.21 
C 

OO(NO)  «  {N0*0.25)  -  0.25 
C 

OSO(NL.NO)  -  8000.  *  EXP(-ALPHA 

c 

900  CONTINUE 
C 

800  CONTINUE 
C 

RETURN 

END 


DD(ND)) 
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